/***************************************************************************************************
		* The Global Persistence of Work from Home

	This code produces Figures 1 and 2 in the paper. 
	
	Figure 1 uses "GSWA_2022to2025.dta" data, with GSWA responses from wave 2 (2022), 3 (2023) and 4 (2024/25) in a panel of 23 countries.

	Figure 2 uses "GSWA_w4" data, with GSWA responses from wave 4 (2024/25) in 40 countries.
 	
***************************************************************************************************/




		**** PATHS:
		
		gl path "C:\Users\pz0511\Dropbox\GlobalWFH\Global WFH 2024-25\PNAS Brief Report\Replication Package/"

		cd "${path}"

		
		
		**** Globals:
		
		* The global below is used in Figure 2. 
		
gl country_order `" `"Australia"' `"Canada"' `"Ireland"' `"New Zealand"' `"UK"' `"USA"' `"Austria"' `"Czech Rep."' `"Denmark"' `"Finland"' `"France"' `"Germany"' `"Greece"' `"Hungary"' `"Italy"' `"Netherlands"' `"Norway"' `"Poland"' `"Portugal"' `"Romania"' `"Spain"' `"Sweden"' `"Türkiye"' `"Argentina"' `"Brazil"' `"Chile"' `"Mexico"' `"Egypt"' `"Nigeria"' `"South Africa"' `"China"' `"India"' `"Japan"' `"Malaysia"' `"Singapore"' `"South Korea"' `"Taiwan"' `"Thailand"' `"The Philippines"' `"Vietnam"'  "'



		gl coefplot_rows_continents  `" c1 = "{bf: Average}"  "'
		local jj = 1
		foreach cntry of global country_order {
			local jj = `jj' + 1
			gl coefplot_rows_continents `"  ${coefplot_rows_continents} c`jj' = "`cntry'"  "'
		}

		
		

********************************************************************************
****** Figure 1: Work from home levels have stabilized since 2023
********************************************************************************

use "Data/GSWA_2022to2025.dta", clear

	* Create region variable:
		* Note that we exclude Brazil from any region
	cap drop region
	gen region = 1 if inlist(country,"Australia","Canada","UK","USA")
	replace region = 2 if inlist(country,"Austria", "France", "Germany", "Greece", "Hungary") | inlist(country,"Italy", "Netherlands", "Poland", "Spain", "Sweden", "Türkiye")
	replace region = 3 if inlist(country,"China", "India", "Japan", "Malaysia", "Singapore", "South Korea", "Taiwan")
	label def region 1 "English Speaking" 2 "Europe" 3 "Asia"
	label val region region



****** Figure: 			
gen all = 1

foreach nn of numlist 1/3 {
matrix values_`nn' = J(9,3,.)

local i = 1
* By overall/gender/region/age_group:
foreach var of varlist all female region age_group {
	qui: levelsof `var', local(levels)
	foreach val of local levels {
		di "`val' in `var'"
		qui: mean n_wfh if wave == (`nn'+1) & `var' == `val'
		matrix values_`nn'[`i',1] = r(table)[1,1]
		matrix values_`nn'[`i',2] = r(table)[5,1]
		matrix values_`nn'[`i',3] = r(table)[6,1]
	
	local i = `i' + 1
	}
}


matrix values_`nn' = (values_`nn')'
}
			
			
			
coefplot (matrix(values_1), mlabgap(3) ci((2 3)) mlabel(string(@b, "%5.2f")) mlabsize(vsmall) se(2) ciopts(recast(rcap)) recast(bar) barwidth(0.25) bcolor(%70)) ///
		 (matrix(values_2), mlabgap(3) ci((2 3)) mlabel(string(@b, "%5.2f")) mlabsize(vsmall) se(2) ciopts(recast(rcap)) recast(bar) barwidth(0.25) bcolor(%70)) ///
		 (matrix(values_3), mlabgap(3) ci((2 3)) mlabel(string(@b, "%5.2f")) mlabsize(vsmall) se(2) ciopts(recast(rcap)) recast(bar) barwidth(0.25) bcolor(%70)) ///
		, coeflabel(r1 = "{bf:Overall}" r2 = "Men" r3 = "Women" r4 = "English Speaking" r5 = "Europe" r6 = "Asia" r7 = "20-33" r8 = "34-46" r9 = "47-64") vertical  legend(order(2 "2022" 4 "2023" 6 "2024/25") pos(6) row(1)) ytitle("Average") ylabel(0(0.5)2) groups(r2 r3 = "{bf:By Gender}" r4 r6 = "{bf:By Region}" r7 r9 = "{bf:By Age Group}", gap(0.5)) xsize(10) ysize(4)

		gr export "Output/Figure1.png", as(png) replace 
		gr export "Output/Figure1.eps", as(eps) replace 

	
	
	* Save the p-values:
matrix Pvalues = J(9,5,.)



* By overall/gender/continent/age_group:
		local i = 1
foreach var of varlist all female region age_group {
	qui: levelsof `var', local(levels)
	foreach val of local levels {
		foreach nn of numlist 1/3 {
		qui: mean n_wfh if wave == (`nn'+1) & `var' == `val'
		matrix Pvalues[`i',`nn'] = r(table)[1,1]	
		}
		* Pvalue between 2 and 3:
		qui: reg n_wfh b2.wave if inlist(wave,2,3) & `var' == `val' , r
		matrix Pvalues[`i',4] = r(table)[4,2]	
		* Pvalue between 3 and 4:
		qui: reg n_wfh b3.wave if inlist(wave,3,4) & `var' == `val' , r
		matrix Pvalues[`i',5] = r(table)[4,2]
			
				di "`var' `val' `i'"

			local i = `i' + 1

	}
}


matlist Pvalues		


frmttable using "Output/Figure1_PValues.rtf", statmat(Pvalues) ///
		ctitle("","(1)", "(2)", "(3)", "(4)", "(5)",  \ ///
			"" , "2022", "2023", "2024/25", " 2023 - 2022", " 2024/25 - 2023" , \  ///
				" ", " ", "", " ", "P-value", "P-value" ) ///
				rtitle("Overall" \ "Men" \ "Women" \ "English Speaking" \ "Europe" \ "Asia" \ "20-33" \ "34-46" \ "47-64" ) replace sdec(3)


				
			
		
************************************************************
**** Figure 2: Working from home levels by country
************************************************************

use "Data/GSWA_w4.dta", clear


mean n_wfh, over(order_countries)


* WFH Levels - College graduates:
	foreach var of varlist n_wfh {

* Mean over all countries
mean `var', over(order_countries)

* We store the coefficients in a matrix called S1_"variable". 
* We will then add the cross-country average, and divide it into many matrices, for graphical purposes
matrix S1_`var' = J(1,40,.)
matrix S1_`var'[1,1.40] = r(table)[1,1..40]

* Now we preserve the data but convert the results to a matrix sequence called results and ress.
matrix results = J(2,41,.)
cap drop regr
cap drop `var'_raw
gen regr = S1_`var'[1,_n] in 1/40

* Add the cross-country average:
sum regr 
ren regr `var'_raw 
gl meanval = r(mean)
matrix results[1,1] = $meanval
matrix results[1,2.41] = S1_`var'[1,1..40]
 

	* Graphical decisions:
 * X-axis title:
if "`var'" == "n_wfh" local ytitles `" "Number of days working from home this week"  "'

* X-axis numeric labels:
if "`var'" == "n_wfh" local ylabels `"0(1)2"'

* Number of the graph:
if "`var'" == "n_wfh" local numm 1

* Number of decimals in the figure:
local decval 1
local mlabel `"  string(@b, "%5.`decval'f") "'



** For graphical purposes, divide results in a matrix in each continent
matrix ress1 = results[1,1] // Average
matrix ress2 = results[1,2..7] // English Speaking
matrix ress3 = results[1,8..24] // Europe
matrix ress4 = results[1,25..28] // Latin America
matrix ress5 = results[1,29..31] // Africa
matrix ress6 = results[1,32..41] // Asia	


* Version 2:
coefplot (matrix(ress1), recast(bar) bcolor(black*0.9)) ///
		 (matrix(ress2), recast(bar) bcolor(cranberry*1.4)) ///
		 (matrix(ress3), recast(bar) bcolor(ebblue*1.4)) ///
		 (matrix(ress4), recast(bar) bcolor(gold*0.6)) ///
		 (matrix(ress5), recast(bar) bcolor(dkorange)) ///
		 (matrix(ress6), recast(bar) bcolor(forest_green*0.75)) ///
		 (matrix(results), mlabel(`mlabel') mlabposition(3) mlabsize(vsmall) mlabcolor(black) mcolor(none) msize(vtiny) msymbol(d) aux(2)),  ///
		legend(off) coeflabels(${coefplot_rows_continents}, labsize(tiny)) ///
		xlabel(`ylabels')  offset(0)  xtitle(`ytitles', size(small)) ///
		ylabel(, labsize(vsmall)) graphregion(color(white)) grid(none) ///
		fintensity(inten90) headings(c2 = "{bf:English Speaking}" c8 ="{bf:Europe}"  c25 = "{bf:Latin America}" c29 = "{bf:Africa}" c32 = "{bf:Asia}", gap(0.3) labsize(vsmall)) 

		* Save graph
		gr export "Output/Figure2.png", replace
		gr export "Output/Figure2.eps", replace as(eps)



		}

				
				
